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Abstract 

The numerical simulation of two-dimensional decaying turbulence in a large but finite 
box presented in this paper uncovered two physically different regimes of enstrophy 
decay. During the initial stage, the enstrophy Q, generated by a random Gaussian 
initial condition, decays as Q(t) oc t~ 7 with 7 ~ 0.7 — 0.8. After that, the flow 
undergoes a transition to a gas or fluid composed of distinct vortices. Simultaneously, 
the magnitude of the decay exponent crosses over to 7 ~ 0.4. An exact relation for 
the total number of vortices, N(t), in terms of the mean circulation T of an individual 
vortex N(t) oc 1/T 2 is derived. A theory predicting N(t) oc and the magnitudes of 
exponents 7 = 2/5 and £ = 4/5 is presented and the possibility of an additional very 
late-time cross-over to 7 = 1/3 and £ = 2/3 is also discussed. 
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The interest in the problem of two-dimensional (2D) turbulence was excited by the pioneering 
work of Onzager 1 , who recognized the role of point- vortices in 2D turbulence dynamics and 
suggested a theory based on some general concepts of statistical physics. In Onzager's 
days the quantitative investigation of two-dimensional hydrodynamics was inhibited by the 
pausity of well-controlled physical experiments. Now, with advent of powerful computers, 
various features of the problem are much better understood due to the remarkable works of 
McWilliams and collaborators 2-5 , Benzi et. al. 6 and many others. Still, many important 
questions remained unanswered. 

If a fluid in a large three-dimensional vessel of linear dimension L is stirred on a length- 
scale £ < i, the non-linear process of generation of the small-scale velocity fluctuations 
leads to energy dissipation, making a steady state possible even in the limit of vanishing 
viscosity. In this case, universality means independence of the small-scale features of the flow 
upon length-scale L. In two dimensions, the situation is very different; the Euler equations 
conserve both energy and enstrophy. As a result, the most prominent dynamic feature of 
2D-turbulence are two fluxes in the wave- number space 7 . The (negative) energy flux leads 
to formation of the slow-varying (low-frequency) large-scale motions containing almost the 
entire energy of the system on a time-dependent integral scale lo(t). The second, enstrophy 
flux, leads to generation of rapidly-varying small-scale velocity fluctuations responsible for 
the enstrophy of the flow. The short-time dynamics, when l (t)/L <C 1, are characterized by 
a quasi-steady inertial range with the energy spectrum E(k) oc k~z and the integral scale 
lo(t) ~ 2ir/ko(t) oc ti <C L. The most interesting phenomenon happens later in the evolution 
when l Q {t) ~ L; the energy spectrum starts piling up at k m n/L forming very slow varying 
and very powerful motions where the major fraction of kinetic energy is concentrated 8-13 . 
The enstrophy associated with these motions is negligibly small. Simultaneously, powerful 
small-scale vortices, containing almost the entire enstrophy of the flow, are formed at the 
centers of these large scale structures. The characteristic scale (radius), a{t) ~ l/kf pa const, 
of these vortices is very small, a(t) <C L. These complex structures depend on geometry of 
the problem and thus are not universal 8 ' 12 . 

The time evolution of decaying two-dimensional turbulence also has three distinct intervals. 
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During the initial stage, the "inverse cascade" leads to population of the scales k < k (t = 0) 
with formation of vortices which merge and create the ever larger vortices. At the same 
time the enstrophy dissipation and some intermittency are generated around small-scales 14 . 
The total kinetic energy in the system remains almost constant while the vortex merger 
persists until only two coherent vortices with positive and negative vorticity are formed 
(McWilliams 2 ). 

Depending on the setup, the experimental and numerical data on the enstrophy decay vary 
rather dramatically. It has been shown that after some transient, the enstrophy Q, in 
turbulence initially generated as a Gaussian random velocity field, decays with time as 
Sl(t) « t~ 7 with 7 « 0.7 - 0.8 15 . The recent theory 16 of 2D-turbulence decay in an infinite 
system led to an asymptotic magnitude 7 = 2/3, which is close to numerical results. Later 
in the evolution, after a substantial depletion of the initial enstrophy, formation of the 
well-separated vortices containing almost entire enstrophy of a flow has been observed 2 . 
This phenomenon inspired various investigators to study the dynamics of 2D-turbulence 
as evolution of a set of initially created point- vortices 3-6 . This set up, however, yielded a 
different decay law for Q(t) with 7 r* 0.35 — 0.4 2-6 . 

To explain this magnitude of the decay exponent, Carnevale et.al. 3 based their scaling theory 
on numerical results by McWilliams 2 giving the total number of vortices in a flow iV r* 
with £ Rj 0.7 — 0.75. More recent simulations by Bracco et al 4 confirmed this result giving 
£ = 0.75 ± 0.05. Observing that the entire kinetic energy, K, of the system is concentrated 
in a rotational motion of the well-separated vortices, the authors of Ref.[3] obtained K rj 
piV(t)c<; 2 a 4 = const and assuming that the peak vorticity of a vortex uj q rj const, they 
obtained their decay law: a(t) oc iV _ 3 Q(t) rj N(t)u>la 2 (t) oc t~^ in a good agreement with 
numerical data. The observed value of the exponent £ ~ 0.75 remained unexplained. The 
theory of Benzi et. al. 6 was based on similar assumptions and led to almost qualitatively 
identical results. Unification of the outcomes obtained from these different setups is one of 
the goals of the work described below. 

In this paper we report the results of a numerical investigation of the 2D turbulence decay 
in a large but finite square. (The definition of "large" is given below.) It is shown that after 
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an initial enstrophy decay, characterized by 7 ~ 0.7 — 0.8, the flow undergoes transition to 
a state that can be represented as a gas of vortices. The enstrophy decay law in this state 
is f2 oc t~ 7 with 7 Rs 0.4. The theory developed in this paper leads to the characteristic 
length-scale of a vortex a(t) oc exactly as in Ref.[3]. The proposed analytic model 

for the vortex interaction gives for the number of vortices in a system as N(t) oc t~ as and 
f2 oc and are in a close agreement with numerical simulations. An exact expression for 
N(t) oc 1/r 2 , where T is a circulation of an individual vortex, is derived directly from the 
Navier-Stokes equations. 

The problem of decay of two-dimensional (2D) turbulence we are interested in this paper 
is formulated in a following way. Consider a time evolution of an initial velocity field, 
u(x, y, 0) = u . defined on a two-dimensional square such that — L < x,y < L. The field u is 
a Gaussian random noise defined in the Fourier space in the vicinity of A; = k (t = 0) = O(^) 
with Ug = 0(1) and iV = const = 0(1). At the initial instant, t — 0, the enstrophy is 



f2 = to 2 = (V x u) 2 = 0(1). The square is large meaning that we are interested in the limit 
v — > and k (t = 0)L — > 00. Still, the box is finite so we will be able to study both short 
time, when k (t)L ^> 1, and long time asymptotics when k (t)L ps 1. In all our simulations 
the initial kinetic energy of the flow K(t = 0) = \p J v 2 (x, t = 0)dx = \ and the fluid density 
p = l. 

The Navier-Stokes equations with the 0(z/ 4 V 4 v) hyper- viscous dissipation terms were sim- 
ulated using a pseudo-spectral method. To be sure that hyper- viscosity did not influence 
the results, two simulations with resolutions 1024 2 and 2048 2 were conducted. The initial 
random field was Gaussian with energy spectrum E(k) = a s -^- (jf) + e ( s+2 K fc p) where 
a * = ^"^1 + } anc ^ s = 3 as in Ref.[15]. The parameters of the simulations are given in the 
following table. 

On Fig.l we show the time evolution of both enstrophy, Q, and energy, K, in a flow. During 
the entire simulation, kinetic energy stays constant; this is not so for the enstrophy. Initially 
after some transient, the enstrophy Q(t) obeys the power law f2 oc t -7 with 7 ps 0.7 — 0.8 and 
is very close to the result reported before by Chasnov 15 . Later in the evolution, when the 
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N 


time-step 


2K 
P 


n(t = o) 


k p 


i/ 4 


1024 


0.0005 


1 


1170 


32 


2.1xl0- 10 


2048 


0.00025 


1 


1170 


32 


9.6xlQ- n 



Table 1: Numerical simulation parameters 



dimensionless time T is such that T = yuJ 2 (t = 0)t > 200, a cross-over to the value 7 ~ 0.4 
is clearly seen. On Figs. 2-3 the time evolution of vorticity field is presented. We can observe 
that the cross-over to a new long-time decay regime (T > 200), with 7 xs 0.4, is accompanied 
by formation of the well-separated distinct vortices. It is interesting, that similar structures 
were observed in case of the small-scale forced 2D-turbulence 9,10 . One feature of Figs. 2-3 
deserves a special attention. At the initial stage of decay, 10 < T < 200, the vortices are 
more or less homogeneously packed, while at the later times the vortex density shows a clear 
tendency of the flow structures to clamp together keeping the mean distance between the 
neighbors more or less constant. This may be an indication of a weakening of the screening 
of the inter-vortex interaction as the vortex density decreases. The merger itself happens 
in two stages; first the same-sign vortices approach each other with translational velocity 
v and start moving as a couple ("vortex molecule") rotating around their mutual center of 
mass. The merging process is relatively slow meaning that during this "pre-merger" period, 
the characteristic frequency of this joint motion is uj m where u m ~ v /a(t) <C u Q and oj is 
the mean angular velocity of an individual vortex. This "dance" is interrupted by a collapse 
with a single emerging vortex with the peak vorticity u Q ~ constant. 



First we present a model which is a generalization of a remarkable work by Carnevale et.all. 3 . 
After that we will demonstrate its consistency with the Navier-Stokes equations. The total 
kinetic energy of a flow ([K] = grew 2 / sec 2 ) is: 

K{t) = P -l v 2 (*)d 2 x = \ P N(t)[v 2 a 2 + u 2 a 4 } = K(t = 0) = 0(1) = constant (1) 
We conjecture and later verify that while the 0(y 2 a 2 ) contribution to the total energy is 
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(a) (1024 2 ) (b) (2048 2 ) 

Figure 1: Enstrophy and energy decay as function of time 

small, it entirely dominates the time-scale of the merger process. According to Carnevale et 
al 3 and our own numerical simulations we assume that uo m const. Neglecting the 0(v%) , 
the equation (1) leads to the estimate obtained in Ref.[3]: \pN (t)uj'l(t)a{t) A = K = const 

and a (0M^fe)*«w-3 

Following Carnevale et al 3 , we introduce the mean distance between the nearest- neighbor 
vortices as R ~ Lj N{t) and estimate translational velocity v Q of a vortex moving in the 
field of other vortices: as: 

v ~ 1 ~^E-^- V (2) 

where V m uo a 2 is a circulation. To complete the argument, we estimate the rotational 
velocity of a vortex as: Vr ~ uj a(t) = ( 2 ^m) ) ^ • Thus, the ratio ^ = (^Tp-)^A^^ oc t~± — > 
and is indeed small. It is important that according to the above estimate, while vr grows 
with time, translational velocity v Q is a time-independent constant. 

Now we need a model for N(t). We consider a vortex merger process as a binary "chemical 
reaction" 2N N with the "reaction rate" TZ(t), where is the inverse of the time the 
vortex pair spends in a bound state prior to the merger event. If vortices are well separated, 
the population equation is: 
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T=1 .72758, min=-1 65.4, max=1 54.3 T=68.436, min=-1 67.8, max=1 54.9 



Figure 2: Vorticity field before crossover (1024 2 ) 




T=205.274, min=-174.6, max=154.4 T=821.043, min=-170.8, max=147.5 



Figure 3: Vorticity field during and after crossover (1024 2 ) 



— = -1lN 2 (t) 
dt y ' 



(3) 



We recall that two same-sign vortices attract each other, tending to form a spinning dipole 
(bound state) with the linear dimensions (radius of the orbit) depending upon the impact 
parameter, relative translational velocity, and other characteristics of the collision process. 
It is important that the larger the relative velocity , the larger the radius of an orbit of a 
bound state. The rotational frequency of this bound pair is 0(v o /I) where / is the impact 
parameter. In the system we are dealing with, the binary collisions happen in the field of 
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Figure 4: Time evolution of energy spectrum (1024 2 ) 

other vortices of both signs often strongly perturbing the pre-merging dipole. It is clear, and 
we observe it in our simulations, that only the most strongly interacting pairs stay stable 
under perturbations from surrounding fluid long enough to eventually merge. The most 
stable ones are those with linear dimension (impact parameter) not larger than a few a(t) 
and having smallest relative velocity. Since the translational velocity v ~ const, describing 
the motion of coupled vortex pair before the merger event, is the smallest velocity-scale in 
the system we choose it as a rate determining parameter. The merger rate is thus estimated 



as: 7Z ~ v D /a(t) ~ v { N ^^ oP )i = ^J^^-N* . Substituting this into the above differential 



2K 

equation gives: 

N = - N ^b'* (4 > 

(£ + l)» <o 

where t ~ \J^t~ an d is an undetermined constant. This gives the enstrophy decay law as 
H(t) « N(t)uila 2 {t) oc t~r, and is extremely close to numerical result. It follows from (4) that 
in an infinite system (L — > oo) the cross-over to 7 ~ 0.4 never happens and the enstrophy de- 
cay process with 7 = 2/3 persists in the limit t — > 00. Now we would like to discuss another 
possibility. If the mean distance between the nearest vortices is large (R = 0(L/ ^jN(t)), 

1 2 1 

then the reaction rate 1Z = 0(v Q VN/L) and equation (5) gives N(t) oc t~3 and £l(t) oc t~s. 
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The relations derived above are the outcome of a kinetic-like considerations. It is interesting 
to verify their consistency with the Navier-Stokes equations. Choosing the displacement 
vector r parallel to the a;-axis and assuming statistical homogeneity and isotropy of a flow 
, the Navier-Stokes equations lead to the relation for the correlation function of the scalar 

• • j— . — 7— r 7 r-16-19 

vorticity = uj{x)uj[x + r) : 

dD w 1 d 

- - = — — —r(v(x + r) — v(x))(ou(x + r) — uj(x)) 2 (5) 
at 4r or 

Multiplying this equation by r and integrating over the interval < r < L gives 16-18 : 

Q r i 

— / drruj{x + r)uj(x) = —-r(v(x + r) — v(x))(uj(x + r) — uj(x)) 2 \q (6) 

In the limit L — > oo the right side of this equation is equal to zero. To evaluate the 
integral, we recall that in 2D hydrodynamics the point-vortex representation of vorticity 
is: uj{x) = X^Tj^x — Xj(t)), where Xj(i) is a position of the i th vortex and Tj stands for 
its circulation. Then, recalling that the flow is filled with vortices of both signs so that 
FiTj = 0, we derive readily: 

J dr m>(xV(x + ri) = J dr^-^^r^ J dx5(x - X;(t))5(x - Xj(t) + ri) = 

V T ^ - V r * - "UP* ~ N ^ a4 ~ canst (7) 
^2^-^2^-^L 2 2^~ C ° nSt (7) 

i,j t=l 

This relation, derived from purely hydrodynamic considerations, is quite important since 
it serves as a correct definition of both oj and a(t) which are the basic ingredients of the 
"kinetic theory" developed in Carnevale et.al 3 . The average- splitting T 2 w uj 2 a 4 , used in 
derivation of the last relation in (7) is an approximation valid only if the peak vorticity and 
vortex radius are statistically independent. At the present time we cannot assess the quality 
of this approximation. 



We can also start with the Navier-Stokes equations and write equation for the correlation 
function of the x-component of velocity field ( here denoted as v ) 16 ^ 18 : 

d 1 d 

— v (x + r)v{x) = 7r^'K~ r ( v ( x + r ) — v(x)) 3 (8) 

C/6 \Ji U I 
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Repeating the considerations leading to (6) we obtain: / r 3 v(x)v(x + ir)dr = const and 
introducing the correlation length a(t), this integral is estimated readily giving N(t)v 2 a 2 ps 
N(t)(v 2 + uj 2 a 4 ) ~ const where iV(i) = L 2 /a 2 . This is exactly the relation (1) evaluated 
above from the energy conservation considerations. 

The main result of this work is unification of the two different exponents 7 ~ 0.7 — 0.8 
and 7 ~ 0.4 as corresponding to intermediate and long time regimes of the enstrophy decay 
in 2D turbulence, respectively. The exact relation N(t) oc 1/T 2 derived from the Navier- 
Stokes equation, is a direct generalization of an approximate formula N(t) oc a -4 obtained 
by Carenevale et.al. 3 using the "kinetic theory" approach and assuming the constancy of the 
peak vorticity of individual vortices. 

We have identified three possible values of exponent 7. According to ref.[16], the initial stage 
of enstrophy decay in a large box (L — > 00) is characterized by 7 = 2/3 which is close to 
results of numerical simulations by Chasnov 15 and the ones presented in this paper. It has 
been demonstrated that in a finite box, after the depletion of the initial externally introduced 
enstrophy, a flow composed of long-living distinct vortices emerges and the exponent crosses 
over to 7 ~ 0.4, very close our theoretical prediction. Another cross-over to 7 ~ 1/3 and 
£ = 2/3 as a very long-time asymptotics in a system with small number of vortices is 
predicted. To the best of our knowledge, this regime has not yet been observed. 
The nature of the transition remains poorly understood. It has been shown 9-10 that in case 
of the small-scale forced 2D turbulence, the well-separated vortices emerge as a result of the 
finite-size effect (condensation). It is hard to say if this conclusion holds for freely decaying 
2D turbulence but if it does, then the question of universality of exponents 7 and £ becomes 
extremely interesting and important. Some indications that the finite-size effects can play 
an important role can be found from Figs. 1-4. It can be shown theoretically that as long as 
ko(t)L ^> 1, the large-scale energy spectrum in this developing flow is E(k <C ko(t)) oc k 3 (fi(t). 
We can see that in the short-time interval, T < 200, this is indeed approximately correct and 
the exponent is 7 ^ 0.7 — 0.8. At the longer times (T > 200), when k (t)L ps 1, the energy 
spectrum starts strongly deviate from the /c 3 -curve. This is definitely due to the finite size- 
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effects. Simultaneously, the well-separated vortices emerge and the exponent crosses over to 
7 ~ 0.4. 
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